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PREDICTION  OF  FUTURE  OBSERVATIONS 
IN  POLYNOMIAL  GROWTH  CURVE  MODELS 
PART  -  1* 
by 

C.  Radhakrishna  Rao 
University  of  Pittsburgh,  Pittsburgh 


ABSTRACT 

The  problem  considered  is  that  of  simultaneous  prediction  of  future 
measurements  on  a  given  number  of  individuals  using  their  past  measurements. 
Assuming  a  polynomial  growth  curve  model,  a  number  of  methods  are  proposed 
and  their  relative  efficiencies  in  terms  of  the  compound  mean  square  pre¬ 
diction  error  (CMSPE)  are  compared.  There  is  a  similarity  between  the 
problem  of  simultaneous  estimation  of  parameters  as  considered  by  Stein 
and  that  of  simultaneous  prediction  of  future  observations.  It  is  found 
that  the  empirical  Bayes  predictor  (EBP)  based  on  the  empirical  Bayes  esti¬ 
mator  (EBE)  of  the  unknown  vector  parameters  in  several  linear  models  pro¬ 
posed  by  the  author  (Rao,  1975)  has  the  best  possible  efficiency  compared 
to  the  others  studied.  The  problem  of  determining  the  appropriate  degree 
of  the  polynomial  growth  curve  is  also  studied  from  the  point  of  view  of 
minimising  the  CMSPE. 


AMS  Classification:  62C12,  62J07. 

Key  Words:  Empirical  Bayes  procedure,  Compound  decision  problem,  Growth 
curves,  James-Stein  estimators,  Ridge  regression. 

*The  paper  is  based  on  a  talk  given  at  the  Indian  Statistical  Institute, 

Calcutta  in  December  1981,  during  the  Golden  Jubilee  celebrations. 
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1.  INTRODUCTION 


Let 
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be  k  linear  models,  where  are  observable  p-vector  random  variables,  6^ 
are  unknown  m-vector  parameters,  X  is  (p x m)  and  x  is  (m*l)  given  matrices. 
The  problem  is  to  predict  y^,...,yk  given  Y  .....Y^  when  the  dispersion 
matrix  of  the  error  term  (£i»r'i)  is  of  the  form 
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are  of  full  rank  and  V ^  f  0.  Suitable  modifications  can  be  made  when  these 
matrices  are  deficient  in  rank  (see  Rao,  1973,  pp.  296-302). 

If  no  assumption  is  made  about  the  joint  distribution  of  (Y^,y^), 
then  the  least  squares  theory  may  be  applied  to  estimate  and  y^  simul¬ 
taneously.  This  leads  to  minimization  of 


(1.3) 


with  respect  to  8^  and  y^.  The  estimates  so  obtained  are  easily  seen  to  be 

Bi  “  (X'V^X)1  X’v'j^ 
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where  8*  is  the  least  squares  estimator  of  6^  from  the  model  i.  The  same  es¬ 
timator  can  be  deduced  by  considering  a  linear  function  L'Y^  +  a  and  min¬ 


imizing  the  mean  square  error 


2 


E(yi‘-’!i"a)2 

subject  to  the  unbiasedness  condition 

E(yi'L’Va)  ’  °*  (1.5) 

Predictors  of  the  type  (1.4),  which  may  be  called  the  best  linear  unbiased 
predictor  (BLUP) ,  have  been  studied  by  Rao  (1973,  p.  234)  and  Toutenburgh 
(1970).  The  MSPE  (mean  square  prediction  error)  of  the  BLUP,  in  (1.4),  is 

E(yryi)2-  °2<Y22~Y2iYn!i2 + 1  (1,6) 

where  d  -  X_x’viiv12  and  u  “  (X'V”*X)-1. 

In  this  paper,  we  examine  the  possibility  of  constructing  predictors 
of  y^  without  using  the  condition  (1.5),  such  that  the  CMSPE  (the  compound 
MSPE) 

Eljj^-y^2]  (1.7) 

is  a  minimum.  Such  a  procedure  leads  to  predictors  analogous  to  Stein  type 
estimators  in  simultaneous  estimation  of  parameters  (see  Stein,  1955  and 
James  and  Stein,  1961) . 


(2.1) 


3 


f  _ 1  W  _  *  It” 


where  d  ■  x-X'V^V^  an<*  M  "  L-V11V12*  T^e  ^irst  tenn  on  c^e  right  hand  side 
of  (2.1)  does  not  Involve  L  and  a,  and  in  order  to  minimize  the  CMSPE  we 
need  only  consider  the  second  term.  It  is  easily  seen  that  the  minimum  is 
attained  for  given  8^,..., 6^  when 


L  -  (Vu+XFX')"1(V12+XFx)  and  a  =  (x'-L'X)ct 


(2.2) 


where 


a  ■  k  ^  I  and  a^F  ■  k  ^KB^-a)  (B^-a) ' . 


(2.3) 


The  optimum  values  of  L  and  a  as  obtained  in  (2.2)  involve  two  functions 

2 

a  and  F  of  the  unknown  parameters  and  a  .  If  a  and  F  are  known, 

then  the  BLP  of  y  which  minimize  the  CMSPE  are 


1-1 . k’ 


(2.4) 
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where  U  *  (X'V^X)  ,  8^  is  the  least  squares  estimator  of  and 


6i  ’  6i  '  uCF+u)~1(6^-a). 


(2.5) 


It  may  be  noted  that  when  a  and  F  are  known,  is  the  best  linear  esti¬ 
mator  of  S^,  i  »  l,...,k,  in  the  sense  that  the  matrix 

k  k  ,  . 

E  ^(^i'?i)(!i‘!i),_E  [(!i"!i)(5i"!i>’ 


(2.6) 


where  bi  are  any  linear  estimators  of  di  and  the  expectations  are  taken  for 
fixed  8^, ...,6^,  is  non-negative  definite.  Further,  may  also  be  recognized 
as  the  Bayes  estimator  of  S^  when  has  the  multivariate  normal  distribution 


4 


and  has  the  prior  multivariate  normal  distribution  with  mean  a  and 

2 

variance-covariance  matrix  a  F,  in  the  sense  that  the  matrix 


ECbj-Bj)  (b^)  )  (bJ-B^ 


(2.7) 


where  the  expectations  are  taken  over  variations  of  Y^  and  8^,  is  non-negative 
definite  (see  Lindley  and  Smith,  1972). 

The  average  MSPE  for  the  BLP's  in  (2.4)  is 

l,- ^  u  r  (A  *  _  ..  \  ^ 


k"  E 


°2(y22-y2ii[Ihi2+f"?) 

-  aVlKF+U)*1!^. 


(2.8) 


Comparing  the  average  MSPE's  of  the  BLUP's  given  in  (1.6)  and  the  BLP's 
given  in  (2.8),  we  find  that  the  last  term  in  (2.8)  represents  the  reduction 
in  loss  when  a  and  F  are  known. 


3.  EMPIRICAL  BAYES  PREDICTORS  (EBP) 

The  best  linear  estimator  6^  of  6^^  as  defined  in  (2.5)  involves  the  know¬ 
ledge  of  o  and  F.  If  they  are  not  known, we  can  substitute  for  them  suitable 

2 

estimates  and  obtain  modified  estimators  of  B^>  Natural  estimators  of  a,  o 
2 

and  G  -  a  (F+U)  are  of  the  form 

A  -1  i 

a  -  k  Z  Bj 


£ 

>  c-1v-c-1  T(Y'v"1Y  -Y'V-1X3i) 
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where  and  c 2  are  constants.  Substituting  these  estimates  in  (2.5)  we 


obtain 


S*  -  B*  -  c  w UB_1(6^-a) 


(3.1) 


where  c  is  a  constant  to  be  suitably  chosen.  The  estimator  6^  which  is  no 
longer  linear  in  Y^  may  be  called  the  EBE  (empirical  Bayes  estimator,  al¬ 
though  the  terminology  may  not  be  appropriate  without  introducing  an  apriori 
distribution  for  8^)  •  Estimators  of  the  type  (3.1)  have  been  considered  by 
Efron  and  Morris  (1972,  1975)  and  Rao  (1953,  1975). 

In  the  paper  (Rao,  1975),  the  author  has  shown  that  when  c  in  (3.1)  is 
chosen  as  (k-m~2)/ (kp~km+2) ,  where  m  *  rank  of  X,  the  expectation  of  the  matrix 


I  (6. -6.) (87-6.) ’  -  l  (B?-B .)  (S?-B . ) ' 

1*3.  -1  '  ~  i*l  -1  -1  -1  -1 

for  any  fixed  B^,...,Bk,  is  non-negative  definite,  which  implies  that 

V.  k. 

ElKd’Bj-d’B.)2]  >  E[I(d’B®-d’B.)2]. 

1  '  ~  '  i  *  ‘  '  -1 


(3.2) 


(3.3) 


g  b 

Now  substituting  8^  for  Bj^  in  (2. A),  we  obtain  the  modified  predictor 


for  y. 


*i  ■  r  ?!  +  Ynhrfi 


which  may  be  called  the  EBP  (empirical  Bayes  predictor) . 


(3.4) 


Upto  now  we  have  not  made  any  distributional  assumptions.  If  we  assume 
a  multivariate  normal  distribution  for  (Y^,y^),  then 


■frjy  -  f'fi  +  YziYuIi 


k  k  k 

Ell^-y^2]  "E  J[y1-E(yi|Y1))2+E  f  [E(y±  | Y±) -y±3 2 

^  ^  ^  k 

-  k  o2(v22-y21v-Jyi2)  +  e  ICd's^-d's)2. 
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The  corresponding  expression  for  Che  BLUP  (1.4)  is 

k  k 

E  I(y1-yi)2-k  °2<Y22"!2i!n!i2)+E  [d’ei"!’?i)2-  (3-6) 

Using  the  result  (3.3),  and  comparing  (3.5)  and  (3.6),  we  find  that 

,2, 


k  .2,  k 


(3.7) 


E[I(y,-yJ  ]  E[£(y,-y.)  j 

1  1  1  1  1 

where  the  expectations  are  taken  for  any  fixed  set  8^ . 8,.  of  the  true 

values.  Thus  the  EBP’s  of  y^....^^  are  uniformly  better  than  the  BLUP's. 


4.  PREDICTION  IN  POLYNOMIAL  GROWTH  CURVE  MODELS 
Let  y  ^  be  the  measurement  at  time  t  on  individual  i.  We  consider  the 
problem  of  predicting  y  ^  ^  on  the  basis  of  y^,...^^  assuming  a  polynomial 
growth  curve  model 


ytl  ■  80i*0<,:)+---+S5i*s<t)+'ti 


t  *  1,...,  p^*  1  j  i  -  1, . . .  ,k 


where  ^(t)  ,i|»^(t)  , . . .  are  orthogonal  polynomials  such  that 


\>  (t)<J»  (t)  ■  1  for  r«m  and  0  for  n^ra 
r  m 


and 


V(€ti)  “  0  *  cov(eti,Cui^  “  °»  t,tu, 
COv(£ti*€uj)  "  °*  i’1^’ 


(4.1) 


Under  the  model  (4.1),  the  least  squares  estimators  of  8^,  o  are 

8^  -  f  yti^r^C^’  i-l,...,k,  (4.2) 

ri  t-i  r 

32  -  ill  *  !l(BrV2l*  k(P“»-D.  (^.3) 
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Also,  for  given  i 

VCe^J-a2  and  cov(0^i , 3^)  =  0 ,  r^q. 

We  consider  different  methods  of  predicting  y  .  . ,  i=l . k,  simultaneously 

p+i  »  i 

and  compare  their  CMSPE  (compound  mean  square  prediction  error) . 

4.1  BLUP  with  a  subset  of  terms 

If  we  choose  only  the  first  (q+1)  terms  in  the  model  (4.1),  then  the 

BLUP  of  y  .is 
p+l,i 

yp+l,i  =  e0i’J/0(p+1)+'--+65i*q(p+1)  ('*.1.D 

and  the  MSPE  for  given  i  is 

°2  it'i'  (p+i)]2  +  [  y  b  ^  (p+i)]2 

0  r  q+1  rl  r 

when  in  fact  the  true  model  has  all  the  (s+1)  terms.  The  corresponding 
CMSPE  is 

2  q  2  k  s  2 
ko  £[*  (p+1)]  +  l  [  I  B  .♦  (P+1)]  .  (4.1.2) 

0  i=l  r=q+l  rl  r 

If  all  the  (s+1)  terms  in  (4.1)  are  used,  then  the  CMSPE  is 

2  s  2 

ko  l  U»r(p+1)]  .  (4.1.3) 

0 

The  omission  of  the  last  s-q  terms  in  (4.1)  provides  better  prediction,  al¬ 
though  the  corresponding  regression  coefficients  may  not  be  zero,  if 
s  k  s 

ko2  l  [*  (p+1)]2  >  l  [l  8  .*(p+l)]2  (4.1.4) 

q+1  r  i-1  q+1  rl  r 

which  might  hold  when  the  last  regression  coefficients  are  small.  The  best 
choice  of  q  is  that  value  for  which  (4.1.2)  is  a  minimum. 
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In  practice,  the  minimization  of  (4.1.2)  over  q  cannot  be  carried  out 
2 

since  o  and  6^^  are  not  known.  An  estimate  of  q  may  be  obtained  by  mini¬ 
mizing  an  estimate  of  (4.1.2),  which  is 

g  k  s  s 

ko2  I[iMp+l)]2  +  l  [  l  6*  *  (p+l)]2-2ka2  £  [*(p+l)]2.  (4.1.5) 

0  i*l  r»q+l  q+1 

Since  the  first  term  in  (4.1.5)  does  not  depend  on  q,  we  need  only  minimize 

the  expression 

k  S  s 

l  t  y  8*  (p+1)]2  -  2ko2  y  [*(p+l)]2  (4.1.6) 

i»l  r**q+l  q+1 

which  is  analogous  to  the  criteria  suggested  by  Akaike  (1973)  and  Shibata 
(1981)  in  the  context  of  fitting  a  model  to  observed  data.  The  emphasis  in 
our  case  is  on  the  prediction  of  future  observations  and  the  procedures  sug¬ 
gested  by  Akaike  and  Shibata  for  obtaining  a  good  fit  to  observed  data  may 
not  be  appropriate. 

Table  1  gives  the  wieghts  of  13  male  mice  measured  at  intervals  of  3 
days  over  the  21  days  from  birth  to  weaning,  as  reported  by  Williams  and 
Izenman  (1981) . 


Table  1.  Weights  of  13  male  mice  measured  at  successive 
intervals  of  3  days  over  21  days  from  birth  to  weaning 


mice' 

^days  3 

6 

9 

12 

15 

18 

21 

1 

0.109 

0.388 

0.621 

0.823 

1.078 

1.132 

1.191 

2 

0.218 

0.393 

0.568 

0.729 

0.839 

0.852 

1.004 

3 

0.211 

0.394 

0.549 

0.700 

0.783 

0.870 

0.925 

4 

0.209 

0.419 

0.645 

0.850 

1.001 

1.026 

1.069 

5 

0.193 

0.362 

0.520 

0.530 

0.641 

0.640* 

0.751 

6 

0.201 

0.361 

0.502 

0.530 

0.657 

0.762 

0.888 

7 

0.202 

0.370 

0.498 

0.650 

0.795 

0.858 

0.910 

8 

0.190 

0.350 

0.510 

0.666 

0.819 

0.879 

0.929 

9 

0.219 

0.399 

0.578 

0.699 

0.709 

0.822 

0.953 

10 

0.255 

0.400 

0.545 

0.690 

0.796 

0.825 

0.836 

11 

0.224 

0.381 

0.577 

0.756 

0.869 

0.929 

0.999 

12 

0.187 

0.329 

0.441 

0.525 

0.589 

0.621 

0.796 

13 

0.471 

0.770 

1.105 

This  could  be  a  recording  error,  but  no  change  was  made 
in  the  present  computations. 


1 

< 


In  each  case,  the  seventh  measurement  is  predicted  using  the  first  six  by  BLUP, 
the  formula  (4.1.1),  for  different  values  of  q  (degree  of  the  polynomial). 

The  sums  of  squared  differences  (SSD)  between  the  observed  and  predicted  over 
the  13  mice  for  each  q  were  as  follows. 

q:  0  1  2  3  4  5 

SSD:  1.7911  0.2063  0.1042  0.1750  0.5991  7.4700 


It  is  interesting  to  note  that  the  second  degree  polynomial  provides  the  best 
model  for  predicting  the  seventh  observation,  and  the  predicion  becomes  worse 
as  we  increase  the  degree  of  the  polynomial,  although  higher  degree  poly¬ 
nomials  should  theoretically  provide  a  better  fit  to  the  data. 

Since  it  is  found  that  a  second  degree  polynomial  is  the  appropriate 
model  for  predicting  future  observations,  we  shall  explore  alternative 
methods  of  estimating  the  regression  coefficients  ^Qi’^li’^i^  for  Predic_ 
tion  purposes  and  examine  their  predictive  efficiencies. 


4.2  James-Stein  regression  predictor  (JSRP) 

III  - 2 

The  least  squares  estimators,  ®oi’^li’®2i  and  °  (with  P  =  b>  s  =  3,  k  =  13) , 
are  computed  as  in  formulae  (4.2)  and  (4.3).  Using  these,  the  J-S  estimators 
of  ®oi’^li’^2i  are  obCa^ned  as  follows: 


,,  26a2. al  c2  l(al  .2 

a - -2-)e  ,  S  -  ICB  )  . 

28S^  0 


(4.2.1) 


The  predictor  of  y^  based  on  the  estimators  (4.2.1)  is 

y7i  "  80iVP+1)  +  euVp+1)  +  S2i*2(p+1)'  (4.2.2) 


It  is  known  that  the  J-S  estimators  have  smaller  compound  mean  square  error 
than  the  least  square  estimators,  which  may  not  imply  that  any  linear  function 
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of  the  unknown  parameters  is  better  estimated  by  substituting  the  J-S 
estimators  instead  of  the  least  square  estimators  (see  Rao  and  Shinozaki, 
1978  and  Rao,  1980). 


4.3  Shrunken  regression  predictor  (SRP) 

_  -1H  i  ~2 

Let  -  k  (Brj+. .  .+8^)  and  a  be  as  computed  in  (4.3)  with  p*>,  s  =  3 

and  k«13.  Then  estimate  Sri  by  shrinking  the  least  squares  estimators 
using  the  formula 


•It 


82 

r 


b2*;2  ri 

r 


. ,  r  =  0,1,2. 


(4.3.1) 


The  predictor  of  y^  based  on  (4.3.1)  is 


y7i  =  60i*0(p+1)+6^1<p+1)+6^°(p+1) 


lirl' 


2i  2 


(4.3.1) 


4.4  Empirical  Bayes  predictor  (EBP) 
Let 


8  -  (6o.8j.6j)'. 


13 


B  -[(sf-8)(B*-B)' 
1  -1  '  ~x  ' 


Then  the  EBE  of  Bi  -  ^eoi’8li’ 62i^  is>  SS  shown  in  Rao  (1975)» 

.  J.  _  (k-q-3) (p-s-l)o2  B-l,Rf_g* 

8i  ?i  k(p-s-l)+2  !  (!i  V' 


(4.4.1) 


The  predictor  of  y^  based  on  (4.4.1),  with  k»13,  q*2,  p«6,  s«3,  is 


y7i  "  B0iVP+1)  +  Sli'('l<P+i)  +  S2i*2(p+1)'  (4.4.2) 

It  may  be  noted  that  8®  is  also  a  Stein  type  estimator  (see  Efron  and  Morris, 

l 

1972  and  Rao,  1975)  of  a  vector  parameter.  Since  the  estimators  6ri»  r* 0,1,2 
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t 


-A 


% 


and  i*l,...,13  are  all  uncorrelated ,  Che  problem  may  also  be  considered  as 
the  simultaneous  estimation  of  the  39  parameters  8ri»  This  was  not  tried  in 
the  present  data  analysis. 

4.5  Ridge  Regression  predictor 

The  ridge  regression  estimator  of  8^  is  computed  from  the  formula  (see 
Hoel  and  Kennard,  1970) 


si 


8 


s2+3a2  si* 

l 


Si  ■  si>2- 


(4.5.1) 


The  predictor  of  y^  based  on  (4.5.1)  is 


7i 


60iVP+1>  +  eii*l(p+1)+62iMp+1) 


(4.3.2) 


The  sums  of  squared  differences  between  observed  and  predicted  values 
for  the  7-th  measurement  over  the  13  mice  for  different  methods  using  a  second 
degree  polynomial  were  as  follows: 


Method:  BLUP  JSRP  SRP  EBP  RRP 

SSD:  .1042  .1044  .0972  ,0951  .1047 


In  the  problem  considered,  ridging  (RRP)  and  James-Stein  procedure  (JSRP)  do 
not  seem  to  improve  the  least  squares  estimators  (BLUP)  for  prediction  purposes. 

As  established  theoretically  in  section  3,  the  EBP  showed  the  best  performance, 
while  SRP  is  a  close  competitor. 

It  is  proposed  to  study  further  methods  of  estimation  of  regression  para¬ 
meters  and  also  transformations  of  the  time  axis  to  improve  predictive  efficiency. 


A 


12 


BIBLIOGRAPHY 


[1]  Akaike,  H.  (1973).  Information  theory  and  an  extension  of  the  maximum 
likelihood  principle.  In  2nd  International  Symposium  on  Information 
Theory,  Eds.  B.N.  Petrov  and  F.  Csaki,  pp.  267-81,  Budapest,  Akademia 
Kiado. 

[2]  Efron,  B.  and  Morris,  C.  (1972).  Empirical  Bayes  on  vector  observations. 
Biometrika  55,  335-47. 

[3]  Efron,  B.  and  Morris,  C.  (1975).  Data  analysis  using  Stein's  estimator 
and  its  generalizations.  J.  Am.  Statist.  Assoc.  70,  311-319. 

[4]  Hoel,  A.E.  And  Kennard,  R.W.  (1970).  Ridge  regression:  Biased  estima¬ 
tion  for  nonorthogonal  problems.  Technometrics ,  12,  55-68. 

[5]  James,  W.  and  Stein,  C.  (1961).  Estimation  with  quadratic  loss.  Proc. 
4th  Berkeley  Symp.  1,  362-379. 

[6]  Lindley,  D.V.  and  Smith,  A.F.M.  (1972).  Bayes  estimates  for  the  linear 
model  (with  discussion).  J.R.  Statist.  Soc.  B  34,  1-41. 

[7]  Rao,  C.  Radhakrishna  (1953).  Discriminant  function  for  genetic  differ¬ 
entiation  and  selection.  Sankhya  12,  229-246. 

[8]  Rao,  C.  Radhakrishna  (1973).  Linear  Statistical  Inference  and  its 
Applications,  Wiley,  New  York. 

[9]  Rao,  C.  Radhakrishna  (1975).  Simultaneous  estimation  of  parameters 
in  different  linear  models  and  applications  to  biometric  problems. 
Biometrics.  31,  545-554. 

[10]  Rao,  C.  Radhakrishna  and  Shinozaki,  N.  (1978).  Precision  of  individual 
estimators  in  simultaneous  estimation  of  parameters.  Biometrika  65, 
23-30. 

[11]  Rao,  C.  Radhakrishna  (1980).  Some  comments  on  the  minimum  mean  square 
error  as  a  criterion  of  estimation.  In  Statistics  and  Related  Topics, 
pp.  123-143  (Eds.  CscSrgB,  Dawson,  Rao  and  Saleh). 

[12]  Shibata,  R.  (1981).  An  optimal  selection  of  regression  variables. 
Biometrika  68,  45-54. 

[13]  Stein,  C.  (1955).  Inadmissibility  of  the  usual  estimator  for  the  mean 
of  a  multivariate  normal  distribution.  Proc.  3rd  Berkeley  Symp.  1, 
197-206. 


13 


[14]  Toutenburgh,  H.  (1970).  Vorhersage  im  allegemeinen  linearen  regression 
modell  mit  stochastishen  regressoren.  Op.  Forshung  und  Math.  Stat.  2, 
105-116. 

[15]  Williams,  J.S.  and  Izenman,  A.J.  (1981).  A  class  of  linear  spectral 
models  and  analysis  for  the  study  of  longitudinal  data.  Tech.  Rept. 
Dept,  of  Stat.,  Colorado  State  University. 


e 


*. 


